setwd("E:/0 公共数据库差异情况/db_for_5hmc/CMC等DEG数据库的分析")
setwd("E:/0 公共数据库差异情况/")
#deg=read.table("sig_DEG.csv",head=T,sep=",")
#degid=unique(deg$symbol)
library(openxlsx)
sci=read.xlsx("2018 Science supp/aat8127_Table_S1a.xlsx",sheet="DGE")
sci=sci[sci$ASD.fdr<0.1|sci$SCZ.fdr<0.1|sci$BD.fdr<0.1,]
sci_scz=sci[sci$SCZ.fdr<0.1,]
sci_bd=sci[sci$BD.fdr<0.1,]

sci_scz=data.frame(symbol=sci_scz$gene_name,SCZ.fdr=sci_scz$SCZ.fdr)
sci_bd=data.frame(symbol=sci_bd$gene_name,BD.fdr=sci_bd$BD.fdr)

scz_0704=read.xlsx("db_for_5hmc/scz-details1.xlsx",sheet = 1)##这是李淼新课题组预测的易感基因
bd_0704=read.xlsx("db_for_5hmc/bip-details1.xlsx",sheet = 1)##这也是李淼新课题组预测的易感基因
scz_0704$minP_IVW=as.numeric(scz_0704$minP_IVW)
scz_0704_sig=scz_0704[scz_0704$minP_IVW<0.05,]
scz_0704_sig=scz_0704_sig[!is.na(scz_0704_sig$Gene),]
scz_0704_sig=data.frame(symbol=scz_0704_sig$Gene,minP_IVW=scz_0704_sig$minP_IVW)
bd_0704$minP_IVW=as.numeric(bd_0704$minP_IVW)
bd_0704_sig=bd_0704[bd_0704$minP_IVW<0.05,]
bd_0704_sig=bd_0704_sig[!is.na(bd_0704_sig$Gene),]
bd_0704_sig=data.frame(symbol=bd_0704_sig$Gene,minP_IVW=bd_0704_sig$minP_IVW)

bd_sig=merge(bd_0704_sig,sci_bd,by = "symbol",all = T)
scz_sig=merge(scz_0704_sig,sci_scz,by="symbol",all=T)
bd_sig_id=unique(bd_sig$symbol)[-c(1:3)]
scz_sig_id=unique(scz_sig$symbol)[-c(1:10)]

write.csv(bd_sig,"BD_sig_deg.csv",quote=F,row.names = F)
write.csv(scz_sig,"SCZ_sig_deg.csv",quote=F,row.names = F)


case=read.table("E:/5hmc_file/2_5hmc_yjp_bam/ASM/bayes_p/bias_AShM_BF_no_motif.txt",head=T,sep="\t")
case=case[case$BF_in_DC>10,]
case=tidyr::separate_rows(case,Gene.refGene,sep=";")
caseid=unique(case$Gene.refGene)
caseid=caseid[!caseid=="NONE"]



case_bd=merge(case,bd_sig,by.x = "Gene.refGene",by.y="symbol")
case_scz=merge(case,scz_sig,by.x = "Gene.refGene",by.y="symbol")
write.csv(case_scz,"SCZ_sig_deg_overlap_5hmc.csv",quote=F,row.names = F)
write.csv(case_bd,"BD_sig_deg_overlap_5hmc.csv",quote=F,row.names = F)



con.file=list.files(path = "E:/5hmc_file/组织特异性表达/",pattern=".csv")
con.fn=con.file[-c(grep("tissue",con.file))]
con.file=paste0("E:/5hmc_file/组织特异性表达/",con.fn)
#############################################################用基因名富集
############################################################SCZ 
a=length(intersect(scz_sig_id,caseid))
b=length(caseid)-a

col_names=c("con.file","case_overlap","case_not","con_overlap","con_not","OR","p.value")
result=data.frame(matrix(NA,1,ncol=7))
names(result)=col_names
for(i in 1:4){
  con=read.table(con.file[i],head=T,sep=",")
  con=tidyr::separate_rows(con,Gene.refGene,sep=";")
  con$unitID=paste(con$Chr,con$Start,sep = ":")
  conid=unique((con$Gene.refGene))
  conid=conid[!conid=="NONE"]
  c=length(intersect(scz_sig_id,conid))
  d=length(unique(conid))-c
  rt_tmp=data.frame(matrix(NA,1,ncol=7))
  names(rt_tmp)=col_names
  rt_tmp[,1]=con.fn[i]
  rt_tmp[,2]=a
  rt_tmp[,3]=b
  rt_tmp[,4]=c
  rt_tmp[,5]=d
  rt_tmp[,6]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$estimate
  rt_tmp[,7]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$p.value
  result=rbind(result,rt_tmp)
}
result=result[-1,]
write.csv(result,"sig_SCZ_DEG_enrichment_by_geneName.csv")

################################################################BD
a=length(intersect(bd_sig_id,caseid))
b=length(caseid)-a

col_names=c("con.file","case_overlap","case_not","con_overlap","con_not","OR","p.value")
result=data.frame(matrix(NA,1,ncol=7))
names(result)=col_names
for(i in 1:4){
  con=read.table(con.file[i],head=T,sep=",")
  con=tidyr::separate_rows(con,Gene.refGene,sep=";")
  con$unitID=paste(con$Chr,con$Start,sep = ":")
  conid=unique((con$Gene.refGene))
  conid=conid[!conid=="NONE"]
  c=length(intersect(bd_sig_id,conid))
  d=length(unique(conid))-c
  rt_tmp=data.frame(matrix(NA,1,ncol=7))
  names(rt_tmp)=col_names
  rt_tmp[,1]=con.fn[i]
  rt_tmp[,2]=a
  rt_tmp[,3]=b
  rt_tmp[,4]=c
  rt_tmp[,5]=d
  rt_tmp[,6]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$estimate
  rt_tmp[,7]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$p.value
  result=rbind(result,rt_tmp)
}
result=result[-1,]
write.csv(result,"sig_BD_DEG_enrichment_by_geneName.csv")
################################################################用unitID富集
a=length(unique(case_scz$unitID))
b=727-a
col_names=c("con.file","case_overlap","case_not","con_overlap","con_not","OR","p.value")
result=data.frame(matrix(NA,1,ncol=7))
names(result)=col_names
for(i in 1:4){
  con=read.table(con.file[i],head=T,sep=",")
  con=tidyr::separate_rows(con,Gene.refGene,sep=";")
  con$unitID=paste(con$Chr,con$Start,sep = ":")
  con1=merge(con,scz_sig,by.x = "Gene.refGene",by.y = "symbol")
c=length(unique(con1$unitID))
d=length(unique(con$unitID))-c  
rt_tmp=data.frame(matrix(NA,1,ncol=7))
names(rt_tmp)=col_names
rt_tmp[,1]=con.fn[i]
rt_tmp[,2]=a
rt_tmp[,3]=b
rt_tmp[,4]=c
rt_tmp[,5]=d
rt_tmp[,6]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$estimate
rt_tmp[,7]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$p.value
result=rbind(result,rt_tmp)
}
result=result[-1,]
write.csv(result,"sig_SCZ_DEG_enrichment_by_unitID.csv")
############################################## BD
a=length(unique(case_bd$unitID))
b=727-a
col_names=c("con.file","case_overlap","case_not","con_overlap","con_not","OR","p.value")
result=data.frame(matrix(NA,1,ncol=7))
names(result)=col_names
for(i in 1:4){
  con=read.table(con.file[i],head=T,sep=",")
  con=tidyr::separate_rows(con,Gene.refGene,sep=";")
  con$unitID=paste(con$Chr,con$Start,sep = ":")
  con1=merge(con,bd_sig,by.x = "Gene.refGene",by.y = "symbol")
  c=length(unique(con1$unitID))
  d=length(unique(con$unitID))-c  
  rt_tmp=data.frame(matrix(NA,1,ncol=7))
  names(rt_tmp)=col_names
  rt_tmp[,1]=con.fn[i]
  rt_tmp[,2]=a
  rt_tmp[,3]=b
  rt_tmp[,4]=c
  rt_tmp[,5]=d
  rt_tmp[,6]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$estimate
  rt_tmp[,7]=fisher.test(matrix(c(a,b,c,d),nrow = 2))$p.value
  result=rbind(result,rt_tmp)
}
result=result[-1,]
write.csv(result,"sig_BD_DEG_enrichment_by_unitID.csv")
